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ABSTRACT 


This  paper  presents  a method  of  selecting  a near-optimal  starting  solution  for 
a large  class  of  discrete  polynomial  approximation  problems  in  the  Lj  norm. 
While  it  is  possible  to  prove  the  optimality  of  these  advanced  starting 
solutions  for  only  a small  class  of  continuous  polynomial  approximation  pro- 
blems, empirical  evidence  indicates  the  starting  bases  will  be  nearly  optimal 
for  a much  larger  class  of  discrete  problems.  This  paper  presents  the  method 
used  to  determine  the  starting  basis  and  a heuristic  justification  backed  by 
empirical  results  supporting  its  use. 


iii 
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I.  INTRODUCTION 

The  investigation  of  a near-optimal  starting  solution  for  discrete  polynomial 
approximation  problems  in  the  norm  resulted  from  a study  to  develop  a 
methodology  for  comparing  and  evaluating  mathematical  programming  software. ^ 

In  this  study  four  approximation  algorithms^ > 3 ,4 ,5  were  tested  on  a large 
set  of  problems  with  diverse  problem  structures  and  characteristics.  Included 
were  80  problems  used  to  test  code  performance  in  the  polynomial  approximation 
of  a continuous  function,  a specific  type  of  approximation  problem.  In  these 
problems,  the  continuous  function  being  approximated,  the  degree  of  the 
approximating  polynomial,  and  the  number  of  observation  points  were 
independently  varied.  In  all  problems,  the  observations  were  equally  spaced 

over  the  same  interval  of  approximation,  [0,1]. 

© « 

When  the  interpolation  points  defining  the  best  solution  were  tabulated, 
with  the  degree  of  the  polynomial  and  the  number  of  observations  fixed,  but 
varying  the  functional  form,  a recurrent  pattern  emerged  which  indicated  that 
the  points  of  interpolation  remained  relatively  stationary  within  the  normal- 
ized interval  regardless  of  the  functional  form  being  approximated.  This 
pattern  suggested  a method  of  selecting  a starting  basis  which  would  be  nearly 
optimal  for  a large  number  of  approximation  problems. 
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After  further  study,  an  analytic  expression  was  identified  that  can  be  used  to 
calculate  a feasible  starting  solution  which,  when  used  in  the  continuous 
analog  of  the  discrete  approximation  problem,  is  the  optimal  solution  in  the 
Lj  norm. ^ Calculating  this  starting  solution  and  using  it  in  the  discrete 
approximation  problem  also  produced  impressive  empirical  results,  both  in 
terms  of  its  relative  distance  to  the  optimal  solution  and  also  in  the  marked 
reduction  in  the  number  of  iterations  needed  to  reach  the  optimal  solution. 
Furthermore,  the  integration  of  this  method  of  selecting  the  starting  basis 
into  existing  algorithms  is  easy.  An  explanation  of  the  success  of  this 
method  can  be  found  in  the  underlying  approximation  theory  and  is  presented  in 
the  next  two  sections.  The  final  section  reports  the  results  of  using  this 
near-optimal  starting  basis  on  a large  set  of  test  problems. 

II.  THE  DISCRETE  POLYNOMIAL  APPROXIMATION  PROBLEM 

The  general  form  of  the  discrete  Lj  approximation  problem  can  be  formulated 
as  follows. 


Given  n sets  of  observations  on  a dependent  variable  y^  and  a single  indepen- 
dent variable  x^ 

{xj^y*}  » 


determine  m+1  parameters 


6 - Bq*  ^1*  B 2 » • • • » ^m  * 


Z 


i-1 


m 4 

7±  ~ l Bj  x,  I 

j-0  1 


which  minimize 
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This  is  the  polynomial  approximation  of  a continuous  function  when  y-j_  = f(x^), 
where  f(x)  is  a continuous  function  defined  over  an  approximation  interval  I. 
Introducing  variables  d and  d , corresponding  to  the  positive  and  negative 
deviation,  the  approximation  problem  can  be  posed  as  the  following  linear 
programming  problem. 


subject  to 


n 

Minimize  £ 
i=l 


+ 

d + d 
i i 


» 


+ raj 

d - d + yi  - [ 8^  x 

i i 1 -i=n  ^ ^ 


0 , 


+ — 

d _>  0 , d^  _>  0 for  i = 1,  2,  3,  . . . ,n 
i 


For  this  class  of  problems  it  is  not  unreasonable  to  assume  the  problems  are 
of  full  rank  (independent  columns)  and  overdetermined  (the  number  of  obser- 
vations strictly  greater  than  m,  the  number  of  3 parameters  in  the  problem). 

A further  condition  necessary  for  one  to  "predict"  the  best  solution  to  any 
problem  of  this  type  is  that  the  solution  be  unique.  When  the  uniqueness  con- 
dition is  satisfied  and  the  number  of  observations  sufficiently  large,  the 
starting  solution  we  suggest  is  very  close  to  the  optimum.  Unfortunately, 
uniqueness  cannot  be  guaranteed  for  certain  approximation  problems. 
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III.  FUNCTIONAL  ANALYSIS  AND  PROPERTIES  OF  Li  APPROXIMATION  OF  A CONTINUOUS 
FUNCTION 

This  section  presents  the  theoretical  basis  behind  a unique  best  L^  approxima- 
tion solution  for  a special  class  of  continuous  functions.  The  section  begins 
by  presenting  several  properties  of  the  polynomial  approximation  of  a continu- 
ous function  in  the  Lj  norm.  These  properties  establish  the  existence, 
uniqueness  and  optimality  of  a solution.  A description  of  a class  of  func- 
tions which  have  these  properties  is  developed  and  it  will  be  shown  that  the 
starting  basis  we  propose  will  be  the  optimal  Lj  solution.  An  alternate 
representation  of  the  optimal  L^  solution,  which  serves  as  a useful  guide  in 
determining  the  "goodness"  of  an  approximation  is  also  presented  as  is  a 
method  for  calculating  the  starting  basis  easily  and  efficiently. 

Since  the  best  Lj  polynomial  approximation  of  degree  m to  a continuous  func- 
tion, f(x),  intersects  the  function  in-  at  least  m+1  points^,  we  will  examine 
the  additional  characteristics  a function  must  possess  to  insure  the  unique- 
ness of  the  optimal  solution.  The  first  assumption  on  the  function  requires 
the  existence  of  a non-trivial  approximation  problem,  i.e.,  f(x)  is  different 
from  the  Lj  approximation  almost  everywhere.  Had  this  not  been  assumed,  any 
m+1  observations  would  have  defined  the  best  Lj  approximation.  By  requiring 
that  f(x)  have  a m+lst  derivative  that  is  continuous  and  non-constant,  a 
nontrivial  approximation  problem  is  guaranteed.  We  characterize  any  such 
function  as  a "higher-order  function"  than  the  polynomial  of  degree  m.  This 
is  closely  analogous  to  the  term  "higher-order  polynomials”  though  more 


general 
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The  uniqueness  of  the  final  solution  and  ultimately  the  optimality  of  the  sug- 
gested starting  basis  is  established  by  making  one  further  requirement  on  the 
derivatives  of  f(x).  If  the  m+lst  derivative,  fm+^(x) , is  nocvanishing  over 
the  approximation  interval,  then  it  will  be  shown  that  the  best  Lj  polynomial 
approximation  to  f(x)  will  intersect  at  exactly  m+1  observations  and,  there- 
fore, will  be  unique.  With  this  condition  satisfied,  one  can  then  find  the 
m+1  observations  which  define  the  optimal  Lj  approximation  to  the  continuous 
problem,  evaluate  the  function  at  these  points  and  enter  these  observations  to 
the  discrete  problem  basis.  Solving  the  corresponding  (m+1)* (m+1)  system  of 
equations  defined  by  this  new  basis  will  locate  a near-optimal  L]_  solution  to 
the  discrete  problem. 

The  following  examples  and  subsequent  theorems  will  solidify  the  theoretical 
underpinnings  bounding  the  number  of  intersections  between  the  continuous 
function  and  the  best  approximating  polynomial.  First  consider  the  problem  of 
locating  all  the  possible  intersections  between  a second-order  polynomial 
P2(x)  and  a first-order  polynomial  p^(x),  (shown  graphically  in  figure  1). 


FIGURE  1 
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Three  possible  situations  exist:  first,  p^(x)  and  P2(x)  fail  to  intersect  in 

[a,  b] ; second,  pj(x)  and  P2(x)  are  tangent  which  will  produce  only  one  point 
of  intersection;  or  finally,  p^(x)  and  P2(x)  intersect  at  two  points.  The 
problem  of  locating  the  points  of  intersection  of  p^(x)  and  P2(x)  is 
equivalent  to  that  of  finding  all  possible  roots  of  a residual  function  de- 
fined as  their  difference  in  [a,  b], 

r(x)  = p2(x)  - PjCx)  = 0,  x e [a,  b]. 

Since  r(x)  is  binomial,  it  can  have  at  most  two  real-valued  roots.  This 
argument  can  easily  be  extended  for  any  two  polynomials,  pm(x)  and  pra+^(x),  of 
degree  m and  m+1  respectively,  with  no  loss  of  generality. 

The  task  of  establishing  the  number  of  intersections  of  pm(x)  and  a polynomial 
of  degree  m+k,  Pni+k^x^ » f°r  k > 1 , or  any  function  of  higher  order  is  slightly 
more  difficult.  Again,  the  behavior  of  the  derivatives  will  determine  the 
number  of  roots  the  residual  function 

r(x)  = Ptn+kCx)  “ Pm(x), 

will  have  over  [a,  b].  By  assuming  a restriction  on  p^k,  it  is  again  pos- 
sible to  limit  the  number  of  roots  r(x)  may  have  to  at  most  m+1  and  therefore 
define  a larger  class  of  problems  which  will  have  a unique  best  Lj  solution. 
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Lemma  1.  Suppose  a function  f(x)  is  continuously  differentiable  at  least  k 
times  in  an  interval  [a,  b].  Suppose  further  that  fk(x)  > 0 (<0)  in  [a,  b]. 

Then  fk”l(x)  is  strictly  monotone  increasing  (decreasing)  in  [a,  b],  and  there 
is  only  one  possible  zero  of  fk~*(x)  in  [a,  b]. 

Together  with  the  following  theorem  it  will  be  possible  to  establish  upper 
bounds  on  the  number  of  real-valued  roots  for  each  successively  lower-ordered 
derivative . 

Theorem  1 [McCormick^].  Suppose  a function  f(x)  is  continuously  differentiable 
at  least  k times  in  an  interval  [a,  b].  Suppose  that  fk(x)  has  q zeros  in  that 
interval  and  that  they  are  known  and  ordered  as  a < zj  < 22  < Z3  ...  <Zq  _<  b. 
Then  fk”*(x)  is  strictly  monotone  in  [a,  zj],  [z^,  zi+ll>  (for  i=l , . . . ,q-i ) , and 
[zq,  b].  There  are  at  most  q+1  zeros  of  fk_1(x)  in  [a,  b],  Specifically,  there 
may  be  one  in  [a,  z^],  one  in  each  of  [z^,  z^+j],  (for  i=l,...,q-l)  and  one  in 
[zq,  b] . If  fk~^(a) *fk~^(zj)  > 0 there  is  no  zero  in  that  interval,  otherwise 
there  is  exactly  one  there.  If  (for  i=»l,  2,...,q-l)  fk_1(z.£)  *fk~1(zi+i)  > 0, 
there  is  no  zero  in  that  interval.  If  fk“^-(Zq)  *fk“^(b)  > 0,  there  is  no  zero  in 
[Zq,  b].  Otherwise  there  is  exactly  one. 

Proof.  Since  fk(z^)  3 fk(z^+j)  * 0 and  there  are  no  zeros  between  z^  and  zi+j 
then  fk(x)  > 0 for  all  x in  ( zit  zi+1),  or  fk(x)  < 0 in  that  interval.  Thus 
the  hypotheses  of  the  previous  lemma  are  satisfied  and  the  appropriate  con- 
clusion follows  for  [z^t  zi+l 1 • The  other  cases  are  identical. 
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Thus,  for  any  function  f(x)  with  a nonvanishing  and  monotonic  (m+l)st  deriva- 
tive Theorem  1 can  be  applied  successively  to  construct  upper  bounds  to  each 
successively  lower-order  derivative  until  a upper  bound  of  m+1  can  be 
established  for  the  number  of  real  valued  roots  to  f(x)  over  the  interval 
[a,b].  The  upper  bound  on  the  number  of  roots  of  the  residual  function  r(x) 
is,  therefore,  dependent  only  on  the  behavior  of  fm+^(x)  since  the  (m+l)st 
derivative  of  the  approximating  polynomial  pm(x)  is  zero  valued  over  the  real 
line , thus 

rm+1(x)  = fm+1(x)  - pm+1(x) 

m 

= fm+1(x). 

Under  the  assumptions  above  and  by  requiring  that  the  associated  approximation 
problem  is  of  full  rank,  the  existence  of  a unique  optimal  solution  to  the 
discrete  polynomial  approximation  problem  defined  by  exactly  m+1  observations 
is  guaranteed.  The  final  question  remains.  How  to  select  the  m+1  obser- 
vations close  to  the  optimal  solution? 

As  shown  in  R±yil  inland  Rice  there  exists  a large  class  of  continuous 

functions  for  which  the  optimal  solution  in  the  continuous  approximation 
problem  can  be  determined  a priori  whenever  it  is  known  a priori  that  the 
optimal  solution  is  unique  and  that  the  interpolation  of  these  two  functions 
occur  at  m+1  points.  The  necessary  and  sufficient  requirements  on  f(x)  for 
this  to  occur  are  1)  f(x)  is  continuous,  2)  the  difference  of  the  function 
f(x)  and  the  Lj  approximating  polynomial  Pm(x)is  different  from  zero  almost 
everywhere  on  [a,b],  3)  the  residual  function,  r(x)-f (x)-Pm(x)  , changes  In 
sign  at  a unique  set  of  m+1  points  in  the  interval  [a,b].  These  points  of 
interpolation  are  determined  in  the  following  theorem. 


9 


Theorem  2:  If  f(x)  is  continuous  and  differentiable,  and  f(x)  and  pm(x) 

intersect  nrt-1  times  in  the  interval  [-1,  1],  then  the  least  approximation 
is  the  unique  polynomial  pm*(x)  which  satisfies 


Pm 


i H ill  , , 

(cos ) = f(cos ),  for  j=l , 2,  m+l.1-1 

m+2  m+2 


til 

This  optimal  set  of  observations  are  the  roots  of  the  mc  degree  Chebyshev 
polynomial  of  the  second  kind,  Um(x) . The  explicit  expression  for  this  family 
of  polynomials  is  listed  below: 


Um(x) 


m/2 

= l (-1 


i=l 


(m-i) ! 
i ! (m-2i ) ! 


(2x) 


m-2i 


or  equivalently 


sin(nH-l)0 

Um(cos  0)  for  x=cos0  ; 

sin© 


which  have  the  following  recurrence  relation 

Urc+iU)  - 2xUm(x)  - Um_ x ( x ) . 1 2 


The  first  four  of  these  polynomials  are: 
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Uo^x)  “ 1» 

UjCx)  = 2x, 

U2(x)  = 4x2-1; 

U3(x)  = 8x3_4x# 

The  series  in  which  these  polynomials  are  used  in  approximation  to  a con- 
tinuous function  f(x), 

f(x)  = lk  bjUj (x) , 
j=0 

are  called  Chebyshev  Series  and  if 

lim  b-j^  = 0 > 

^ ->-00 

the  series  is  termed  a Chebyshev  series  expansion  of  the  function  f(x). 

Chebyshev  series  expansions  have  long  been  known  to  converge  very  quickly  to 
the  target  function  compared  to  other  series  approximations  and  can  be  a use- 
ful tool  in  determining  the  significance  of  the  neglected  term  in  the  poly- 
nomial approximation.  Generating  the  equivalent  Chebyshev  representation  to 
the  best  Li  approximating  polynomial  provides  the  user  with  information  on  the 
relative  significance  of  the  last  term  of  the  Chebyshev  representation.  When 
the  coefficient  of  the  higher  order  terms  are  very  close  to  zero,  the  user  can 
assume  little  accuracy  will  be  gained  by  increasing  the  degree  of  the  approxi- 
mating polynomial.  Conversely,  if  the  last  term  of  the  Chebyshev  approxima- 
tion are  of  the  same  order  as  the  previous  coefficients,  attempts  to  Include 
more  terms  may  provide  a better  approximation,  or  give  the  user  an  Indication 
that  the  target  function  f(x)  may  not  have  a converging  Chebyshev  series 
expansion  and  an  alternate  approximation  method  should  be  used. 
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IV.  GENERATING  THE  STARTING  BASIS 

The  method  used  in  calculating  the  starting  basis  involves  determining  the 
zero-valued  points  of  the  Chebyshev  polynomial.  Using  the  trigonometric 
representation  of  the  Chebyshev  polynomials, 

Sin(m+1 ) 9 

Um(x)  " Um(cos9)  = 

Sin9 

it  is  clear  that  the  zero  values  of  Um(x)  occur  at 

kir 

x.*  = cos  for  k = 1,  2,  . . . , m and  -1  < x^  < 1 

m-t-1 

In  our  work,  the  zero-valued  points  of  Um(x)  were  translated  into  the  interval 
of  approximation,  evaluated  by  the  continuous  function,  and  entered  as  con- 
straints to  the  problem.  Results  of  using  these  points  as  a starting  basis 
are  presented  in  the  next  section. 


V.  COMPUTATIONAL  RESULTS 

In  our  experiment,  256  approximation  problems  were  used  in  testing  the  perfor- 
mance of  the  advanced  starting  method  described.  All  problems  were  produced 
by  an  Li  polynomial  approximation  test-problem  generator1^  using  eight  differ- 
ent continuous  functions,  varying  both  in  form  and  difficulty  (see  Table  la  of 
the  appendix).  The  problems  generated  ranged  in  degree  from  3 to  10,  with 
100,  200,  400,  and  800  observations  over  the  closed  interval  [0,1],  Each 
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problem  was  solved  twice.  One  run  of  all  256  problems  did  not  use  the 
advanced  solution  method  but  included  together  with  the  sorted  equidistant 
observations,  the  m+1  generated  observations  {xj_  } used  in  the  new  starting 
solution  method,  thereby  guaranteeing  identical  problems  in  both  runs.  The 
256  problems  were  then  rerun.  The  observations  (xi  } now  formed  the  first  m+1 
constraints  to  the  problem  and  were  selected  to  enter  the  starting  basis  by 
the  Li  code  on  the  first  m+1  pivots.  The  method  used  to  solve  all  problems 
was  double-precision  L^-approximation  algorithm  developed  by  N.  N. 

Abdelmalek. ^ 

In  this  section,  the  results  from  solving  the  problem  with  the  new  starting 
solution  method  are  compared  to  results  from  the  original  method.  Two  perfor- 
mance measures  were  used  in  the  comparison.  One  measured  the  distance  be- 
tween the  objective  function  value  Z at  the  near-optimal  starting  point,  and 
the  optimum  objective  function  value  Z , normalized  by  the  distance  from  the 

g — * 

starting  objective  function  value,  Z (where  3=0),  to  Z . Or,  mathematically, 
the  normalized  distance  is 

Zi  “ 

D . = i=l,2,  ...,256 

Zs  - Zi* 

A second  measure  was  difference  in  iteration  count  between  using  the  near- 
optimal  "starting  basis"  and  the  starting  basis  defined  by  3=0. 
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The  results  in  measuring  the  performance  of  the  advanced  solution  method  by 
the  normalized  objective  function  distance  are  given  in  Table  2 of  the  appen- 
dix and  summarized 
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FIGURE  2.  Relative  Distance  from  the  Optimum 


Of  the  total  256  problems,  only  12  were  more  than  1 percent  from  the  optimum. 
Closer  examination  of  the  outliers  reveals  that  all  outliers  of  magnitude 
greater  than  1 percent  are  related  to  problems  produced  by  two  of  the  eight 
continuous  functions  being  approximated.  Closer  investigation  reveals  that 
the  distance  D is  worst  in  the  approximations  made  by  polynomials  of  lower  de- 
gree. The  functions  (illustrated  in  Figures  3 and  4)  connected  with  these 
outlier  values  are  highly  nonlinear  and  approximations  by  lesser  degree  poly- 
nomials are  very  poor. 


in  Figure  2. 

Advanced  Starting  Solution  Results 


*The  regular  pattern  of  the  outliers  in  Figure  2 illustrate  that  problems 
of  a specific  function  and  degree  are  the  same  relative  distance  from  the 
optimum  for  100,  200,  400,  and  800  observations. 


Observed  Values 


Y = EXP((6*X  + 1)*SIN(6*X  + 1)) 


Y = (1 50 * X)/(1  + ( 1 5 * X) * *2) 


0.0  0.2  0.4  0.6  0.8  1.0 


Interval  of  Approximation 

FIGURE  3.  FUNCTION  3: 


Interval  of  Approximation 


FIGURE  4.  FUNCTION  5: 


Though  selection  of  a low  degree  polynomial  was  unwise  by  the  user,  the  start- 
ing basis  did  reduce  the  number  of  iterations  required  to  solve  these  problems 
in  13  of  the  16  cases,  as  shown  in  Table  1.  Similar  data  for  problems  of  200, 
400,  and  800  observations  are  given  in  Table  2a  of  the  appendix. 


DEGREE  OF  10  9876543 
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FUNCTION 

#3 

28 

14 

42 

27 

19 

14 

24 
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23 

24 

19 

13 

11 

14 

15 

FUNCTION 

#5 

26 

15 

19 

14 

24 

8 

18 

7 

14 

6 

19 

17 

13 

10 

23 

31 

TABLE  1.  ITERATIONS  NEEDED  TO  OBTAIN  AN  OPTIMAL  SOLUTION 
(NUMBER  OF  OBSERVATIONS  =>  100). 
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Increasing  the  degree  of  the  approximating  polynomial,  in  these  outlier  cases, 
did  reduce  the  values  substantially,  with  4 of  the  16  problems  listed  above 
having  an  advanced  starting  objective  function  value  identical  to  the  optimal, 
to  8 significant  digits. 

Over  all  256  problems,  a total  of  62  problems  using  the  advanced  method  were 

— 8 

at  the  optimum,  requiring  no  additional  pivots  for  an  accuracy  of  10  . A 

total  of  5558  iterations  were  necessary  to  solve  all  256  problems  using  the 
standard  method.  Solving  the  256  problems  using  the  starting  method  reduced 
the  total  number  of  iterations  needed  to  solve  all  problems  by  2028,  a reduc- 
tion of  36  percent.  In  all,  241  problems  had  a reduced  iteration  count,  2 
marked  no  change,  and  13  increased. 

In  problems  where  the  degree  of  the  polynomial  was  overspecified,  a degenerate 
problem  much  more  difficult  to  solve,  is  produced.  We  noted  several  problems 
of  this  type  which  had  an  increased  iteration  count  but  also  had  a very  small 
normalized  distance  measure  This  suggests  there  are  many  other  solutions 

with  objective  function  values  very  close  to  the  optimal  solution.  This  could 
force  the  problem  solver  to  test  a large  number  of  solutions  of  nearly 
identical  objective  function  value  before  obtaining  the  optimal  objective 
function  value.  Solving  these  problems  with  a polynomial  of  lesser  degree 
improves  iteration  count,  and  in  several  cases  required  no  additional  pivots. 
From  this  we  can  infer  that  in  the  problems  where  a lesser  degree  polynomial 
fit  the  function  quite  well,  the  Chebyshev  expansion  converged  very  quickly, 
thereby  producing  a good  approximation  to  the  problem.  Inferences  drawn  from 
the  results  such  as  those  described  above  may  provide  the  user  insight  into 
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the  problem  concerning  "goodness  of  fit”  and  "sufficiency  of  degree"  of  the 
approximating  polynomial. 

The  results  of  this  computational  study  indicate  ways  in  which  L} -approximation 
codes  used  for  polynomial  approximation  problems  can  be  improved.  They  are 
collected  below. 

1.  Use  the  zeroes  of  the  Chebyshev  polynomials  as  the  observations  in 
the  starting  basis. 

2.  Convert  the  Li  polynomial  coefficients  into  the  coefficients  of  the 
equivalent  Chebyshev  polynomial  representation  to  provide  the  user  with  infor- 
mation on  the  rate  of  convergence  of  the  Chebyshev  terms  and  the  correctness 

of  the  user's  degree  specification. 
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8.  Y = CBRT(x) 
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